Simulating DNLS models 



Mario Mulansky 

University of Potsdam, Department of Physics and Astronomy 
Louisiana State University, Center for Computation and Technology (mmulanskyacct .Isu.eduj 

m 

O 

Abstract 

5-H 

M^ We present different techniques to numerically solve the equations of motion for the 

"^ widely studied Discrete Nonlinear Schrodinger equation (DNLS). Being a Hamiltonian 

l/^ system, the DNLS requires symplectic routines for an efficient numerical treatment. Here, 

we introduce different such schemes in detail and compare their performance and accuracy 

by extensive numerical simulations. 
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O 1- Introduction 
O 

5/2 The Discrete Nonlinear Schrodinger Equation (DNLS) has been known to physicists, 

O biologists, chemists and mathematicians for more than 30 years now. It describes a simple 

c/j model of coupled unharmonic oscillators defined by the complex valued phase/amplitude 

^~> tjjn G C at lattice site n. The dynamics of the oscillators is governed by the DNLS that 
writes in one spatial dimension: 



d 

, i-^lpn =Vn1pn+4>n+l+i'n-l+ Pl^nfi^n, (1) 

QO where i is the imaginary unit, Vn G M is the local potential and /? G M the nonlinear 

^^ strength. These equations of motion can be derived from the following Hamilton function: 



/3,,,. ,4 



.• 7j = ^KiV'„r+V':+i^n+V'«+iC + |iV'«r- (2) 

o 

^^ The Hamiltonian character immediately gives an integral of motion which is usually 

. . called energy and denoted as E := H = const. However, the DNLS possesses another 

conserved quantity, the norm J\f := '^^ I'/'nP, usually set to unity N = 1. The existence 
of two conserved quantities makes this equation particularly challenging for numerical 
5_^ approaches, as will be explained later. 

Cd This model first appeared as a description of polarons by Holstein [5T]. Later, the 

DNLS was used by Davydov in his studies of protein dynamics |T3l HI] as well as in 
the context of local modes in small molecules |48| . In recent years, however, two new 
applications of the DNLS gained increasing attention: Bose-Einstein-Condensates (BEC) 
and coupled optical wave guides. 
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In the context of Bose-Einstein condensates, the DNLS describes the mean field ap- 
proximation of a weakly interacting Bose gas at zero temperature as shown e.g. in |35j . 
The predictions of this description have been, to some extent, experimentally verified 
for condensates in optical traps and a periodic potential, e.g. OlS]. A very interesting 
phenomenon that can be studied in BECs is Anderson localization, which denotes the 
"absence of diffusion" [3] in disordered systems. Several experimental realizations of this 
effect have been reported, e.g. [1711151 H?! HI |35] . The corresponding DNLS accompanied 
by a random local potential is often called the Discrete Anderson Nonlinear Schrodinger 
Equation (DANSE) and has been subject of heavy numerical studies in the past years. 

Moreover, the DNLS gives a very accurate description of the propagation of light in 
optical waveguides with a nonlinear material [25l|TT]. As for BECs, optical waveguides 
were also used to study Anderson localization by imposing a disordered potential. Espe- 
cially the fact that the nonlinear strength can be controlled quite precisely in terms of 
the properties of the nonlinear material makes this system very appealing for studying 
Anderson localization and nonlinearity experimentally. This has been done recently by 
Schwartz et al. in [47j . A more extensive review on the history and applications of the 
DNLS can be found in [M]. 

The increasing experimental accessibility of systems described by the DNLS has led to 
numerous numerical investigations of such models. Special attention has been given to the 
case where Anderson localization, induced by disorder, is accompanied by nonlinearity. In 
terms of the DNLS this can be modeled by choosing the local potential Vn to be a random 
variable, typically as independent and identically distributed (iid) from some interval of 
size W: Vn G [—W/2, W/2]. A large number of numerical investigations studied the slow 
destruction of Anderson localization by the nonlinear interactions between the localized 
eigenmodes [131 [TH [Ml UHl [Ml [51] . In all these works, subdiffusive spreading of initially 
localized modes is reported up to computationally accessible times, recently confirmed 
by an experimental study |32j . Additionally, similar spreading laws have been observed 
in quasi-periodic nonlinear lattices [5S] and the nonlinear Stark ladder [T71 [53] . Recent 
results on the scaling properties of chaos [H] and on the scaled spreading in a reduced 
model 3S] raised some questions on the asymptotic validity of the observed subdiffusive 
spreading. Several attempts to find a theoretical description of the spreading have been 
proposed, mainly based on an effective noise theory pjjj, but a full understanding of 
the interplay between disorder and nonlinearity is still lacking [15j. Lately, the theory 
of chaotic diffusion has been developed [37] and successfully applied to a different class 
of systems with disorder and nonlinearity [101 [S] , but can possibly be extended to the 
DNLS/DANSE models as well. Here, however, we are not going to further address 
the very interesting problem of asymptotic spreading, but rather focus on the different 
algorithmic techniques that are used to obtain those numerical results. 

This article is organized as follows: After this introduction, we will review the general 
properties of symplectic integrators and the idea of operator splitting to construct such 
methods in section [2] In section [3] we introduce the different numerical approaches for 
simulating the DNLS model and in section [4] we compare their performance. Finally, we 
end with our conclusions in section [s] 



2. Symplectic Schemes 

Solving the DNLS model here means to find an approximate solution 4'{t) oi the 
equations of motion M from an initial condition tp{t — 0). Being a Hamiltonian system, 
the DNLS has a symplectic flow map and exhibits conservation of energy E — const. 
This requires the usage of symplectic routines that preserve the symplectic nature of 
the system (SH HU] . Formally such a solution can be written in terms of the Liouville 
operator e*^^, defined via Poisson brackets [18] and acting on the initial condition 'ip{t) — 
e^^^iplO). If the action of this operator is known explicitly the system is called integrable 
and the solution 4'{t) can be written analytically. In most cases, however, this solution 
can only be approximated by numerical methods. 

A very common technique to find a symplectic scheme for a non-integrable system is 
the so-called operator splitting [34]. For this, the Hamilton function has to be separable, 
which means it can be written a,s H = A + B, where the action of the operators e*^" and 
e*'^^ are known explicitly. Then, the integration of the solution from time t to t + At for 
some small time step At can be approximated by: 

gAtL„ ^ gAt(L^ + LB) ^ ^AtL^^AtL^ ^ ©(Ai^). (3) 

The splitting scheme above is the most simple and only accurate in first order of Ai. In 
a more general way, the approximation can be written as a product of j operators: 



gAtLij ^ 



■Qga.AtL^gb.AtLB_^0(^^P+l)_ (4) 



The operators e"''^*^-* and e^''^*^^ represent the exact integrations according to Hamil- 
tonians A and B over times at At and biAt respectively. The parameters ai, hi are chosen 
such that the resulting product is an exact representation of e^*'"" of order Ai^. Such a 
splitting scheme effectively computes a trajectory not of the original system with Hamil- 
tonian H = A + B, but of a new Hamiltonian K — A + B -]- ©(At^). Thus, the symplec- 
tic property of the dynamics is preserved. Moreover, the computed trajectory, obtained 
from many subsequent steps Q, conserves a new energy E := K = E + 0{AtP), which 
means that along this trajectory the original energy is not exactly conserved, but only 
fluctuates around its exact value with a magnitude of O(AiP). With a non-symplectic in- 
tegration scheme this is not the case, as there the numerical error of the energy increases 
at every step by an amount 0{At^), which accumulates along the trajectory. Hence, 
especially for long integration times the symplectic routines are essential as they allow 
reasonable energy conservation even for large time steps Ai which makes them much 
more efhcient than non-symplectic methods. Thus, symplectic integrators have become 
standard for Hamiltonian systems and several splitting routines of different orders have 
been developed in recent years, see e.g. [S2 [3 113 UHl El EO] ■ 

In the following, we will use two second order (p = 2) symmetric splittings: the 
SBABi and SBAB2 scheme [30]. Symmetric splittings are especially interesting as they 
naturally lead to self-adjoint algorithms. The time evolution according to those splittings 
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anction H = A + B are defined as follows: 




^AtLH , 


. e^*/2LBeAtL^gAt/2LB +0(Ai3) 


(5) 


K -^ 


= A + B + OiAf-\B\) 




^AtL„ _ 


_ ^6iAtLB^aiAtL^g62AtLi3gaiAtLAgfciAtLi3 _|_ 0C/\^3) 


(6) 



K ^ A + B + 0{At^ ■ \B\ 



with oi = 1/2 , &i = 1/6 and 62 = 2/3 in the last line and K denoting the formal 
Haniiltonian integrated exactly by the given scheme. The symmetric composition ensures 
that these schemes are identical to their inverse with negative timestep: SBABi/2(^i) = 
(SBABi/2(— Ai))"-*^, a property denoted as self-adjoint. 

3. Numerical Methods 

Here, we will present different algorithms to numerically calculate approximate tra- 
jectories of the DNLS equation. At first we focus on methods based on the two-part 
operator splitting described above, all of which have a different approach to solve the 
non-local term in (n]). The first method uses Fourier transform for the non-local part and 
has been applied to the DNLS in (SH HSl EH |24] , it will be abbreviated as "FT" here. 
The second algorithm employs an implicit Crank-Nicolson scheme and has been used 
in [3S1 [32] to integrate the DNLS, it will be referred to as "CN" throughout this text. 
The third method is the so-called "PQ" scheme introduced in ^5 , where the non-local 
part is again separated into two integrable operators. Finally, a new multi-symplectic 
scheme is introduced based on the ideas of Bridges and Reich [6] that was already applied 
to a DNLS without local potential [54] and will be called Euler-Box "EB" scheme here. 

3.1. Fourier Method 

The obvious way to treat the coupling term in the DNLS is to use a spectral method. 
This involves forward and backward Fourier transforms at every time step and is thus 
called "FT" in this text. It has already been applied to the DNLS as described in [^ . 
For this scheme, the Hamilton function (pi) is split into two parts H = Aft + BpT with: 

Aft = ^ V'n-iC + V'«C-i 

n 

Here, A contains only the coupling and B consists of the linear and nonlinear local 
potential. This implies that the action of the time evolution operator Q^^^Apj, becomes 
local when applied to the Fourier transform of the state tjjq. Specifically, e^^^^pj, g^(,^g 
as follow [5T] : 

N 

n=l 

^, ^g-2^Atcos(2.(9-l)/JV)^^ (8) 

n=l 
4 



Q^t^Ap^ 



where the first and the last operation are forward- and backward Fourier transforms of 
the state tpn and the second step is the time evolution in Fourier space. The action of 
^AtLBpr^ can be written explicitly as it is fully local in real space: 

Having found two explicit representations of the two parts of the time evolution one 
can now construct integrators by consecutively applying e -^^t and q^^^^ft , The 
simplest form is given by the first order approximation e^^^" = e^*^^^Te^*^^j^T + 
0{At^). As described above, symmetric second order schemes can be obtained by using 
the SBABi or SBAB2 splitting (Is]), (|6|. A numerical study on the behavior of the energy 
error for the simple first order and the two second order splitting schemes of the FT 
method is shown in Figurefl] We compared the three splitting methods by integrating the 
DNLS system (fTl) with a random potential chosen independent and identically distributed 
Vn G [—2, 2], a nonlinear strength oi /S — 1 and N — 128 lattice sites - the DNLS with a 
random potential is usually referred to as DANSE model. As initial condition we chose 
a Gaussian in the lattice center with width a = 10, ensuring JV = Yl IV-'nP = 1- We 
performed the numerical time evolution using each of the above schemes up to the time 
T = 100 and repeated for decreasing step sizes At = 1.0 .. . 10~^. As quantification for 
the accuracy we calculate the mean square energy error AE and norm ZW along the 
trajectory: 



AE = ^{{E„, - EoV) = /^ J2(E„, - ^o)"' AM=^/{{Kn-Mor) (10) 

y rn 

where the index m denotes the time and Eq is the exact energy and A/q = 1 the exact 
norm, fixed to unity in these simulations. Figure IT] shows the behavior of the energy error 
for the different splittings for decreasing stepsize. One clearly observes the second order 
behavior of the SBABi and SBAB2, compared to the first order trivial splitting ^. 
Additionally, the SBAB2 splitting shows a significantly smaller energy error than the 
simpler ABC and SBABi schemes. For stepsizes smaller than At ^ 10~* the finite double 
precision starts limiting the accuracy of the computation, hence smaller stepsizes At < 
10"'* are not advisable for the SBAB1/2 FT schemes. Note, that the FT based splitting 
schemes exhibit numerically exact norm conservation so no analysis of AJ\f is required. 

3.2. Crank-Nicolson (CN) 

As the Schrodinger equation is of such great importance, many numerical methods 
have been developed to find trajectories in this system. One of the most important ap- 
proaches is the "Crank-Nicolson" scheme jT2] , a second order symplectic scheme to treat 
the usual (linear) Schrodinger equation, also used for solving other partial differential 
equations like the wave equation |44| . To use this scheme, we will split the Hamilton 
function (pi) into a linear and a nonlinear part and then use the ideas of Crank and 
Nicolson to solve the linear part. Hence, we write H = ^cn + ^cn with: 
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Figure 1: Behavior of the energy error for the FT scheme and different sphtting methods. 
The simulation was done for a DNLS with N = 128 lattice sites, /3 = 1 and a random 
potential Vn G [—2,2]. AE is the averaged energy error along the trajectory (10) up to 
a total integration time T — 100. 



The Liouville operators for the two individual terms are then defined by their actions 
on '0 as: 






Vnlpn) 



(12) 

(13) 



The time evolution operator e ^^cn is nonlinear, but can then be written explicitly due 
to its local character: 



e'^*-^BcN 



: i'n 



-iAt^|i/i„|' 



tpn 



(14) 



where ip' denotes the time evolution of one time step At according to the Hamilton 
function B only. 

Unfortunately, c^'^-^cn can only be written explicitly for a constant or periodic po- 
tential. As we want to treat the generic case here, we will use the Crank-Nicolson scheme, 
a second order, norm preserving approximative method. Applying it to the linear part 



La (12 1 (we skip the index CN her for simplicity), the exponential is approximated by 
Caley's formula: e'^*^^ sa jzri~Atl2' ''^hich leads to an implicit equation for the evolved 
state ip': 



^AtL/ 



{l^LAAt/2)i}' = {l + LAAt/2)iP. (15) 

From (12 1 can be seen that La is a band matrix with ~iVn on the diagonal and — i 



on the upper/lower sub-diagonal, this means ( 15 ) is a set of linear equations which 
can be solved by a GauB-algorithm. Moreover, due to the band structure of La the 
computational requirement is only linear in the system size. 

Using the symplectic integrators for the two parts ^cn and -Bcn one can now con- 
struct a symplectic scheme as described in section [21 Namely, we have implemented the 
usual Crank-Nicolson splitting scheme of order p = 1: CN after (pi, the second order 
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Figure 2: Behavior of the energy error for the CN scheme and different spHtting methods. 
The simulation setup is as described in Figure IT] The dashed hne shows exemplarily the 
second order behavior and is located at the same position as in Figure IT] for comparison. 

CNgBABi scheme following (l5| and the CNgBAB2 scheme from ^k. Remember that all 
these schemes are symplectic and norm preserving. We also note that for weak nonlin- 
earities the last method CNsbab2 represents a considerable improvement over CNsbabi, 
as the error gets significantly smaller when |-Bcn| ^1- In a typical spreading setup as 
chosen in the numerical examples of this article, the total error in the energy reduces by 
about one order of magnitude, as seen in Figure [2J 

The results are presented in Figure [2] in terms of the energy error A_B, as the norm 
again is preserved exactly. The plots illustrate the respective p = 1,2 behavior of the 
methods, as well as the improvement of the CNgBAB2 over the more simple CNsbabi- 
For both methods we found an optimal stepsize At w 10^^ that leads to minimal energy 
errors of IS.E w 10~^^. For stepsizes below these values, the finite double precision again 
leads to a linear increase of the error. Furthermore, the error of the CNgBAB2 is clearly 
smaller than for the FTsbab2 shown in Figure [T] as easily seen from comparing with the 
dashed black lines that are plotted at the same position in both graphs. 



3.3. PQ-Splitting 

While in the methods above Fourier transforms or the Crank-Nicolson scheme were 
used as an integrator of the linear coupling term in H, Bodyfelt et al. presented a new 
idea to treat the nonlocal term - the PQ scheme [5]. They proposed to apply another 
operator splitting of the linear part such that the final evolution operator consists of 
three separate mappings. 

To derive this method, we first have to divide the state into real and imaginary part: 
i^n ^ o-n + ibn • The Hamilton function ([2]) is then: 



H^Y^ Vn{al + hi) + 2(a„a„+i + hnh 



'n+l 



) + f(«^ 



bif 



(16) 



This can then be split into the three integrable parts: 

n 



Q ^y ^ 2a„a„+i (17) 






P=2j26„6„+i. 

n 

As above, the time evolution operator can now be approximated: 

PQ : c^^tH = ^AtLs^AtLQ^AtL^ ^ 0{At^), (18) 

where with defining q;„ = 2Vn + (3{a^^ + 6^) the individual steps are given as: 

AtLg . \(^'n = o,n cos(a„ Ai) + 6„ sin(a„ Ai) 
[6^ = 6„ cos(a„At) - a„ sin(a„Ai) 

^AtLQ.l< = an ^20) 

[K^bn~ 2Ai(a„_i + a„+i) 

gAtip . ia'„ = a„ + 2At(6„_i + 6„+i) 

This concatenation of three solvable evolution operators defines a first order symplec- 
tic scheme which is an integrator of the formal Hamiltonian K = B + P + Q + 0{At). 
However, as pointed out in section [21 the error of the splitting can be improved by using 



better approximations than (18). The simplest symmetric second order splitting is given 
by: 

PQabc : e^*^" = e^^-e^^-e^*^«e^^-e^^- + 0{At') (22) 

K = B + P + Q + OiAf ■ \B\\P\\Q\). 

Surprisingly, this obvious three-operator concatenation has been discussed only very 
recently by Skokos et al. [50] where it also was applied to the DNLS. It was called 
"ABC" scheme there, so we denote it as PQabc here. Besides that, one can also use the 
two-operator splitting presented above and apply that two times. This idea was used 
in [S] where the Hamiltonian H is split in two steps as follows: first H = A + B followed 
hy A = P + Q. Applying the Leapfrog method SBABi to the two splittings we get the 
following scheme: 

PQsBABl : e'^*^" = e^Ln^^Lp^^LB^AtLQ^^LB^^Lp^^LB _,. 0(Ai3) (23) 

K = B + P + Q + OiAf ■ \B\\P\\Q\). 
Applying the SBAB2 method individually to the two splittings also results in a second 
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Figure 3: Behavior of the norm error for different splitting methods for the PQ scheme 
and the Euler-Box method. The simulation setup is as described in Figure[TJ The dashed 
line shows exemplarily the second order behavior and is located at the same position as 
in Figures 111 and [2| for comparison. 



order symmetric scheme given by the following 13 successive simple mappings: 

PQ . gAt_H _ gdiAtLi5gdiC2AtLpgC^AtLQgd2C2AtLpgC^AtLQ 

X grfiC2AtLpgd2AtLi3gdiC2AtLj3gC^AtLQ 
^ d2C2AtLp clAtLq diC2AtLp ^diAtLB , Q(^f^) 

K = B + P + Q + 0(Ai2 . iBplPplQI) 



(24) 



with di = 1/6, d-z = 2/3 and C2 = 1/2 [30 . This is again symplectic and has order p = 2. 



For small B this would lead to reduction of the error. However, for the splitting (17), the 
norms of the operators B , P , Q are all ^ 1. Hence, it is not obvious that this scheme 
leads to a better error behavior, but a numerical investigation shows indeed that the 
energy error AE decreases by more than one order of magnitude when using the SBAB2 
method compared to SBABi. This is shown in Figure [31 where we again integrated the 
DNLS (fTl) as above with N = 128 sites, an idd. random potential Vn = —2 ... 2 and 
/3 = 1 up to the time T = 100 using decreasing time steps At — 0.5 . . . 10~^. As initial 
conditions we again used a Gaussian at the center with width a — 10. Because the 
norm is not exactly conserved by the PQ scheme, we will here use AJ\f as quantification 
for the accuracy of the method for the sake of variety. One finds that the PQabc ^^'^ 
the PQsBABi produce very similar errors. The Euler-Box scheme, presented in the next 
chapter, exhibits a slightly smaller error, but by far the best splitting is the PQsbab2- 
Comparing Figures [2] and [3] one finds that the CN based schemes still produce the smallest 
errors. However, the PQ schemes are computationally much simpler which allows smaller 
time steps At at the same CPU time compared to the CN methods. Hence, a detailed 
performance study is required to identify which method is more efhcient. This will be 
done in section ID 

9 



3.4- Multi-symplectic Euler-Box Scheme 

The above schemes rehed on the sphtting of the Hamihon function into integrable 
parts. In this section, we will present a different approach that is based on a multi- 
symplectic formulation of the equations. Multi-symplectic schemes were introduced by 
Bridges and Reich [6 for Hamiltonian PDEs and they are known to locally conserve 
symplectivity. As the DNLS system can be thought of as a discretized version of a 
continuous nonlinear Schrodinger equation, the application of the multi-symplectic theory 
is reasonable [23j. Actually, this has been done already by Wang et al. [51] where a 
multi-symplectic scheme for the DNLS without local potential (T4 = 0) was developed. 
Here, we will adopt this scheme closely following the calculations in [S^ for the general 
DNLS (IT]). As result, we will obtain a locally implicit, self-adjoint multi-symplectic 
scheme of second order. 

To derive this scheme we have to start at the continuous wave equation for the complex 
valued, space- and time-dependent wave function il^{x,t): 



ii't =- i^xx + 27A + V{x)4, + /^IV'I V- 



(25) 



Note, that from a straight forward spatial discretization with grid size Ax ~ 1 one 
immediately obtains (fTl). By setting ip{x,t) = a + ib and v — a^,, w = b^ this can be 
formulated as: 

Mzt + Kz^ = \',S{z), 

with 

^— (a +b) + -{v +w) + -{a 



S = 



by 



(26) 



(27) 



and M and K being anti-symmetric matrices and z containing the state: 
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This formulation emphasizes the multi-symplectic nature of the system. By introducing 
a spatial and timely discretization Xn,n — 1,2, . . . and tm, rn = 1,2, .. . with grid size Ax 
and time step At the discretized state is defined z™ :— z(x„,tm)- The multi-symplectic 
Euler-Box scheme for this discrete state variable is then given as: 



Af+J+z™ + Af_r z™ + i^+5+z™ + K^S^C - V.5(z™), 



(29) 



S^ are the forward/backward difference operators acting on time and space, 

z™)/Ai and (5_z™ = (z™ - z™_i)/Ax. M+, Af_ and K+, K_ are 



m /^?n+l 



i^: 



with 6' 

e.g. (S+z 

splittings of the symplectic structure matrices M — M^ + M_ and K = K^ + K_ , where 

the conservation of symplectivity demands Mj_ ~ — Af_ and K"^ — — A'_. This matrix 

splitting is not unique. However, in [51] it was found that all possible splittings lead to 

one of two fundamental schemes. The first one is given by taking Af+ and K^ as upper 

triangular matrices: 



Af + = 



/O -1 0\ 




, K+ = 


/O -1 




0\ 

-1 




\0 Oyi 




1^0 


o; 



(30) 
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Submitting this into (29) we obtain, after substituting v^ and w™, the fohowing scheme: 



C+' = K ^(C+i - 2a:r + <-i) - Ai(K + 2)a™ 



?7i\2 I /Lrn\2\ ^m 



M^{{a::Y + (C)')a: 



Ai 



C+' = < + ^^(C+\' - 2C+' + C-^-i') + At(K + 2)a:r+' 



(31) 



This is a first order, multi-symplectic scheme [Ml [22]. The scheme is locally implicit 
as the second equation for a™+^ is given implicitly and requires to solve a quadratic 
equation. However, no implicit dependence on neighboring lattice sites is included and 
hence the scheme does not involve solving a set of nonlinear equations which makes it 
locally implicit. From here on, we will set Ax = 1 which makes this method solving the 
DNLS ^. 

The second fundamental splitting is obtained by exchanging M+ o M_ while leaving 
K^ and K_ as is. This gives another scheme similar to the one above but with a locally 
implicit equation for fe™+^ instead of a™^^. It can be shown that the two schemes are 
adjoint to each other with respect to time reversal [M]. The adjoint method of a one-step 
method z'"+-'^ = <I>At^'" is defined as the inverse original map with reversed time —At, 
that is $J^j :— <&lXt and thus can be obtained by solving the following equations for 
^m+i. ^m _ (j)_^j(2;™+-'^). Given the two adjoint schemes above by composition theory 
a new second-order, self-adjoint method can be created by constructing the mapping 
^*At/2^^t/2 [SB, which is given by the following equations: 



m+7 

an ■ 



= <- ^{K+1 - -^K + K-i) At(K + 2)6™ 



- Atp{{a^:Y + {h^Y)h-: 

6r ' - C - ^(C+i' - 2ar ^ + C+Z) - Ai(K + 2)ar ^ 

-Atp{{an^'^f + {bTh^)aT'^ (32) 

a™+i = a;"+^ - ^(C++i' - 2C+' + C+i') - At(K + 2)6™+i 

- Atp{{a^+'f + {V:+'f)U:+\ 

This scheme consists of two explicit (first and third line) and two locally implicit steps 
(second and fourth line) and defines the multi-symplectic Euler-Box (EB) scheme. The 
behavior of the energy error for this scheme is plotted in Figure Is] (black stars) . As ex- 
pected, it is second order accurate, and the error lies between the PQsbabi and PQsbab2 
schemes. However, the Euler-Box method requires to evaluate two root functions for each 
lattice site at each step, which is computationally quite demanding. So again, only a de- 
tailed performance study can reveal if this method presents some advantage over the 
other schemes introduced above. It should be mentioned that this scheme is still a par- 
ticularly interesting approach for solving the DNLS, as it can be easily generalized to two- 
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(a) Fourth order Yoshida composition. 



10' 



10^ 



Ai 



EB6 
CN6 

-At" 



10" 



(b) Sixth order Yoshida composition. 



Figure 4: Energy error of the fourth (a) and sixth (b) order composition methods based 
on the SBAB2 Crank-Nicolson (CN4/CN6) and the Euler-Box (EB4/EB6) scheme. The 
numerical setup is the same as for the results above. 



or even higher-dimensional lattices. A clear advantage over the other methods described 
above, for which a generalization to higher dimensions is not straight forward. 

3.5. High- Order Schemes 

All the methods above were accurate to at most second order. Of course, one often 
requires higher order methods, especially when a good accuracy is wanted. Fortunately, 
there exists a generic way of constructing higher order methods from symmetric, second 
order schemes such as the SBABl, SBAB2 or the Euler-Box algorithm [3T]. Many such 
schemes have been developed, for example by Yoshida [53] , but here we restrict ourselves 
to presenting two higher-order methods, one of order four and one of order six. Assume 
we have a symmetric second order scheme defined by the mapping ^' — ^(At)tjj with 
the stepsize At as parameter. In this context, a scheme is called symmetric if it is self- 
adjoint $(Ai) — <i>*(Ai), where the adjoint scheme is the inverse mapping with negative 
stepsize: $*(At) = <i>^^(— At). For the SBABl/2 schemes it can be easily seen that they 
are indeed symmetric, but this also holds for the Euler-Box scheme, by construction, as 



explained in section 3.4 and in 



Following [S3], we construct a fourth order scheme $^(Ai) as follows: 



$'*(Ai) = 

with ao = -21/3/(2 - 2^/3) and ai 
can be constructed: 



$(Q:iAt)$(aoAi)$(aiAi), (33) 

— 1/(2 — 2^/'^) 53J. Similarly, a sixth order scheme 



$^(Ai) = $(/i3Ai)$(^2Ai)$(//iA<)$(/ioAi)$(^iAt)$(^2Ai)$(^3Ai). 



(34) 



The parameters /io...3 are not defined uniquely in this case. Here we will use the set 
called "solution A" in [531 : 



^il = -1.17767998417887 
Ai3 = 0.784513610477560 



H2 = 0.235573213359357 
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(a) Performance of second order methods. (b) Performance of higher order methods. 

Figure 5: Performance of the different second order schemes (a) and higher order schemes 
based on CNsbab2 (b). These graph show the CPU time to integrate a system with 
N — 1024 lattice sites up to T = 10^ with a given averaged energy error AE on an Intel 
Core 17, 2.93GHz machine. In (b), additionally the results for a non-symplectic 8-th 
order Runge-Kutta scheme (RK) are plotted. 



Note, that these methods of constructing higher order schemes based on the PQabc 
scheme have recently been discussed by Skokos et al. in [50], where they also used the 
DNLS with random potential as benchmarking example. In Figure[4]we examplarily show 
the energy error for the schemes $"* (a) and <1>^ (b) based on the symplectic CNgBAB2 
and the Euler Box schemes. The 0{At'^) and 0{At^) behavior is clearly visible, again 
with a smaller error for the CN based schemes. 



4. Performance 

After having introduced and discussed several numerical methods to solve the Dis- 
crete Nonlinear Schrodinger equation, we will now study the efhciency of the different 
approaches. We will first consider only the second order schemes, as the most efficient of 
these can then also be used to construct the most efficient higher order scheme. For all 
of the methods described above the computational effort increases linearly with system 
sizflM hence no study on different system sizes is required to compare the different meth- 
ods. Note especially that although the Crank-Nicolson (CN) scheme involves solving a 
linear system of N equations, which formally is of complexity 0{N'^), one can make use 
of the band structure of the coupling term Lacn with only three non-zero diagonals and 
thus implement this algorithm with complexity 0{N) as well. 

We analyze the performance of the different methods by computing a trajectory 
for a lattice with N = 1024 sites and a random potential Vn G [—2, 2] and nonlinear 
strength /? = 1. This is the typical setup for studying spreading in the DANSE model [4^ 
1161 I39j . We again start the time evolution from an initial Gaussian with width a = 



^The FT method has an NlogN behavior, so it might become more disadvantageous for very large 
systems 
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10 and total norm A^ = 1. We compute a trajectory up to a time T = 10^ with 
decreasing stepsizes At. As before, we use the mean squared energy error AE as given 



in ( 10 ) as quantification of the accuracy, but here the energy is only computed every 10 



timesteps, so the major computational effort lies in the time evolution. The results for 



the different second order schemes are shown in Figure 5a where the CPU time is plotted 
in dependence of the mean energy error AE. The FTsbab2 and the EB exhibit a similar 
performance, while the PQsbab2 scheme introduced by Skokos et al. [51j represents a 
clear improvement. However, one gets the best performance when using the CNgBAB2 
scheme. It requires about the same computational effort as the FTgBAB2 scheme, but 
its different splitting gives a much smaller energy error, as already seen from comparing 
Figures fl] and [2] We also checked the performance of other splittings (SBABi and ABC), 
but the SBAB2 versions were always superior. 

In Figure [5b] we compare the performance of higher order schemes based on the most 
efficient second order scheme CNsbab2- Namely, we show results for the Yoshida 4 and 



Yoshida 6 compositions introduced in (331 and (34 1. For comparison, we also plot the 
second order results for CNsbab2- Interestingly, the higher order schemes are only more 
efficient when a high accuracy is required. If an accuracy AE > 10~^ is sufficient, the 
second order scheme CNgBAB2 provides the best performance. If one requires smaller 
errors, a higher order scheme should be chosen. In Figure |5b| also the performance of a 
non-symplectic 8-th order Runge-Kutta scheme with stepsize control [20] is presented, 
as implemented in the Boost. Odeint library [1]. We found that for the total integration 
time chosen here, T — 10^, this non-symplectic scheme is competitive and even the most 
efficient for very high accuracy. However, one has to keep in mind that this scheme is 
non-symplectic, hence the energy error AE will increase linearly in time. Thus, for very 
long time scales T > 10^, the symplectic schemes will outperform the controlled Runge- 
Kutta scheme. For the DANSE model that was used here examplarily, the integration 
times go up to T = 10^ i36^ , where the higher order symplectic schemes will surely beat 
the RK method. 

5. Conclusions 

We introduced and described several numerical schemes to compute approximate 
trajectories for the Discrete Nonlinear Schrodinger equation ([l]). The DNLS is a very 
important model with many applications, hence identifying the most efficient numerical 
scheme is of great interest. For our numerical performance test we relied on the DANSE 
setup that is used extensively in the past to study the interplay between disorder and 
nonlinearity. Of the mostly used second order splitting schemes, the CNsbab2 method 
was found to show the best performance in terms of the least CPU time for a given energy 
error AE. Hence, we conclude that this scheme is superior to the others and should be the 
first choice for numerically treating DNLS models when a moderate precision is required. 
A new and particularly interesting approach is the Euler-Box scheme as it can also be 
used for higher dimensional systems. But the implicit treatment of the nonlinearity in 
this case strongly restricts the possibility to treat different local nonlinear potentials (e.g. 
higher powers). The FT, CN and PQ splitting schemes, on the other hand, can be used 
for any local nonlinear potential in the DNLS and are thus more flexible in this regard. 
Note, that all presented schemes here are suitable to additionally integrate the set of 
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linearized equations to obtain Lyapunov exponents. For the CN scheme this has been 
aheady done in [36,. 

Finally, we showed that when a high accuracy is required, higher order schemes 
provide better performance than the standard second order schemes. At an integration 
time of T sa 10^, even non-symplectic schemes are competitive. But this changes with 
increasing T. We note that performance results are always to be taken with care as they 
might vary greatly for different compilers, operating systems or hardware. However, we 
believe that our results are a helpful guide for which simulation code to choose when 
dealing with DNLS models. 
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